Time Domain Mapping of Spin Torque Oscillator Effective Energy 



Graham E. Rowlands, 1 Jordan A. Katine, 2 Juergen Langer, 3 Jian Zhu, 1 and Ilya N. Krivorotov 1 

1 Department of Physics and Astronomy, University of California, Irvine, California 92697, USA 
^Hitachi Global Storage Technologies, 3403 Yerba Buena Road, San Jose, CA 95135 
Singulus Technologies, 63796 Kahl am Main, Germany 
(Dated: February 20, 2013) 

Stochastic dynamics of spin torque oscillators (STOs) can be described in terms of magnetization 
drift and diffusion over a current-dependent effective energy surface given by the Fokker-Planck 
equation. Here we present a method that directly probes this effective energy surface via time- 
resolved measurements of the microwave voltage generated by a STO. We show that the effective 
energy approach provides a simple recipe for predicting spectral line widths and line shapes near the 
generation threshold. Our time domain technique also accurately measures the field-like component 
of spin torque in a wide range of the voltage bias values. 

PACS numbers: 75.70.Cn, 75.75.-c, 75.78.-n 



Spin torque (ST) from a direct spin-polarized current 
[iHj] can excite magnetization auto-oscillations in the 
free layers of nanoscale spin valves and magnetic tunnel 
junctions and thereby generate microwave power 

[ill Il2j at a frequency tunable by the current (l3 - 16 1 



Such spin torque oscillator (STO) devices show potential 
for applications as tunable nanoscale microwave sources 
and magnetic field sensors for computer hard drives 17 1 . 
Due to the STOs' nanoscale dimensions, their auto- 
oscillatory magnetization dynamics are strongly affected 
by thermal fluctuations [Lsl [lij ] , and quantitative under- 
standing of these stochastic dynamics is crucial for the 
development of devices with desired properties such as 
narrow generation line width and high frequency agility. 

In this Letter, we demonstrate a method of using time- 
domain measurements of an STO's voltage oscillations 
(2(i 21 1 for quantitative description of the underlying 
stochastic dynamics. In contrast to frequency domain 
techniques that probe the dynamics indirectly via mea- 
surements of the STO spectral properties, our method 
offers a direct look at time evolution of the magnetiza- 
tion vector. We measure the statistical ensembles of the 
angles at which the magnetization trajectories cross the 
sample plane, and compare them to predictions made 
by theories of stochastic magnetization dynamics. Our 
work demonstrates that the ST-dependent effective en- 
ergy Fokker-Planck formalism [22, 23 [ gives an accurate 
description of the observed ensembles. Based on this ef- 
fective energy approach, we develop a simple recipe for 
predicting spectral line widths and line shapes near the 
generation threshold. This technique also allows us to 
accurately measure the field-like component of ST (FLT) 
over a wide range of voltage biases. 

While the methods discussed here are expected 
to be general, we focus our present study on STOs 
based on magnetic tunnel junctions (MTJ) patterned 
into 150x70 nm 2 elliptical nanopillars from a Ta(5)/ 
PtMn(15)/Co 7 oFe 30 (2.3)/Ru(0.85)/Co4oFe4oB 2 o(2.4)/ 
MgO(0.85)/Co 2 oFe 6 oB 2 o(1.7)/Ta(5) multilayer (thick- 
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FIG. 1: (Color Online) (a) Schematic of the MTJ nanopillar 
with the layer thicknesses given in nm. (b) Measured bias 
dependence of the MTJ conductance in the parallel (Gp) and 
antiparallel (Gap) states. Also shown are the average conduc- 
tance Go and the inverse time-averaged resistance {R}^ 1 for 
the indicated external in-plane field, (c) Examples of trajec- 
tories followed by the free layer magnetization on the sphere, 
and the distributions of the sample plane crossing angles at 
a non-zero temperature. We use spherical coordinates 9 and 
tp, where 8 is the polar angle defined with respect to the sam- 
ple plane normal z and <p is the azimuthal angle defined with 
respect to the P state direction x. (d) A time-domain re- 
sistance trace, R(t), which maxima and minima are used to 
reconstruct the sample plane crossing distributions. 



nesses in nm). Prior to patterning, the multilayer is 
annealed for 2 hours at 300° C in a 1 Tesla in-plane 
magnetic field that sets the pinned layer exchange bias 
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direction parallel to the long axis of the nanopillar. The 
free layer in these structures, pictured in Fig. [lja), 
posseses a large perpendicular magnetic anisotropy 
(PMA) energy E± = K x sin 2 6 + K 2 sin 4 6 that reduces 
both the critical current I c and the frequency of the 
STO auto-oscillations 12J, |24|. Additionally, the free 



layer exhibits voltage controlled magnetic anisotropy 
(VCMA) (ijj 26 1 . The corresponding anisotropy field 



is Hx z = [H p o + AH pa V + H p \ sin z 9\ cos 9, where 
H p o is the first order anisotropy field and AH p q is the 
coefficient of the VCMA field linear in voltage bias. We 
include H p x, the second order anisotropy field, since 
it becomes important due to partial cancellation of 
the out-of-plane shape anisotropy and the first order 
PMA 27]. The combined effect of PMA and VCMA 
is, nevertheless, insufficient to overcome the easy-plane 
magnetic shape anisotropy and the easy axis of the free 
layer magnetization remains in the sample plane. 

In order to extract information on the free layer mag- 
netization trajectories from the STO voltage signal, we 
find the time-varying total resistance across the MTJ, 
R(t), which is written as the sum of time-averaged (R) 
and time-dependent AR(t) components. These compo- 
nents are read out, respectively, by a DC voltmeter and 
a 12 GHz bandwidth 40 GS/s oscilloscope connected to 
the appropriate ports of a bias-T [28[ . To ensure that the 
RF signal far exceeds the noise floor of the oscilloscope 
(5.2 mV rms ), it is amplified by a 35 dB amplifier with a 
low noise figure of 1.3 dB . We assume that the angular 
dependence of the conductance across the MTJ is [29J 
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G = G a (V) + -AG(V)m-p 



(1) 



where G Q (V) = (G A p (In- 
ductance and AG{V) — 



-Gp(V))/2 is the average con- 
G P (V) ~ G AP {V) is the full 
scale conductance change. Here Gap(U) and Gp{V) are, 
respectively, voltage-dependent conductances in the an- 
tiparallel and parallel states of the MTJ, while m and 
p are, respectively, unit vectors in the direction of the 
fixed and free layers' average magnetizations. We obtain 
G Q {V) and AG(V), which are plotted in Fig. EQV), by 
extracting resistance extrema from hysteresis loops of re- 
sistance versus magnetic field taken in the vicinity of the 
easy-axis as described in the supplemental material fijj ]. 
The time-dependent component of the resistance is 

v(t) son + R cx + (R) 



AR(t) 



(2) 



/ 5on 

where / is the DC current applied across the device, 
R cx is extrinsic resistance contribution from contacts and 
sample leads, and V(t) is the voltage signal measured at 
the 50f2 oscilloscope with the microwave circuit amplifi- 
cation and attenuation factored out. 

Since the MTJ resistance depends only on the projec- 
tion of m onto the polarization vector p, we cannot recon- 
struct three-dimensional magnetization trajectories (or- 
bits) from the electrical signals. The orbits are, however, 



symmetric about the sample plane, as seen in Fig. [He). 
Therefore, those points at which m crosses the equator 
correspond to extrema in V(t) (and hence AR(t)) as pic- 
tured in Fig. [XTd). Since the polar angle 8 = n/2 is 
known at these crossings, one may determine by means 
of Eqs. ([TJ and ^ the azimuthal crossing angles 



AG(Vi) \R 



G (Vi) 



(3) 



Here R4 are the extremal resistance values and Vi = IRi 
are the corresponding voltages across the sample. These 
plane crossing angles are histogrammed separately for 
crossings at maxima and minima, and are plotted in Fig. 
[2] for several values of /. These distributions give new 
insights into the STO dynamics. As / is raised, the sep- 
aration between distributions of maxima and minima in- 
crease as the ST increases the amplitude of magnetization 
precession. The individual distributions broaden and be- 
come asymmetric due to changes in the ST-dependent 
effective magnetic energy. The crossing point of the left 
and right distributions (corresponding to the energy min- 
imum at in-plane angle ipo) is observed to shift towards 
the AP state with increasing current. This effect arises 
exclusively from the FLT, as there are no other voltage 
dependent fields acting in the plane of the sample. 

In order to validate the results of our mapping pro- 
cedure, we confirm that the w distributions are repro- 
duced by spin torque theory [J Q ■ The simplest means 
of generating the expected angular distributions is rote 
integration of the stochastic Landau-Lifshitz (LL) equa- 
tion in the macrospin approximation. We make use of a 
graphics processing unit (GPU) to carry out these calcu- 
lations for various realizations of the thermal field. The 
speedup afforded by this method (a factor of at least 10 2 
compared to simulations performed with CPU) allows us 
to fit the macrospin results directly to the experimen- 
tal distributions [28[. The strengths of the in-plane ST 
(which pulls m in the direction m x (p x m)), FLT (an 
effective field along — p), and H p q are taken to be fitting 
parameters. We include H p q as a global fitting parame- 
ter since it exhibits a significant sample-to-sample varia- 
tion presumably arising from free layer inhomogeneities. 
The demagnetization tensor N is assumed to be that of 
an elliptical disk 0, |3(|. All additional input param- 
eters, including the Gilbert damping parameter a and 
the magnitude of the VCMA field, AH p q, are taken from 
independent measurements [261 ] . 

The probability distributions of the plane crossing an- 
gles and the generated microwave signal power spectra 
obtained from these simulations are shown in Figs. [2] 
and[3]Ja). The simulations corroborate our assumption 
that the current-dependent shift of the distributions is 
uniquely determined by the FLT, and we are thus able to 
robustly extract its bias dependence as seen in the inset 
of Fig. [2J In lieu of the fitting method mentioned above, 
we may calculate the FLT directly from the experimental 
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FIG. 2: (Color Online) Plane-crossing probability distribu- 
tions Pc(<p) measured at several values of the bias current for 
H = 500 Oe applied at ipH = 120°. Also shown are distri- 
butions calculated from macrospin and micromagnetic simu- 
lations as well as by the Fokker-Planck approach described 
in the text. The macrospin distributions have been shifted 
towards cp = by 0.03-7T as a visual aid. This shift is an ar- 
tifact of the macrospin approximation [2|| . The inset shows 
the magnitude of the FLT extracted (using both macrospin 
fitting and the Stoner- Wohlfarth model described in the text) 
from the current-induced shifts of the equilibrium angles. 



data in Fig. [5] by using a simple Stoner- Wohlfarth (SW) 
model [Hj]. The FLT is thereby 

Hm = -H x + H y cot(<ys + Sip) 

+4irM s (N X x~ N yy )cos(<po+6tp), (4) 

where H x and H y are the in-plane components of the ex- 
ternal field (including the dipolar contribution) , iV^ and 
N yy are the in-plane components of the demagnetization 
tensor, ipo is the equilibrium angle in the absence of FLT, 
and Sip is the deviation caused explicitly by FLT. Thus, 



by recording the deflections of the crossing distributions 
(at their intersection points), we find that, as shown in 
the inset of Fig. [3J the FLT exhibits an approximately 
quadratic bias dependence of similar strength to that re- 
ported in Ref. [31. In contrast to ST ferromagnetic res- 
onance methods 32|, 33 1, for which careful subtraction 
of background signals is needed at non-zero current bias 
31], Eq. Q provides a fast and simple way of measuring 
FLT strength in a wide range of voltage bias values. 

We perform micromagnetic simulations of the STO dy- 
namics in order to confirm that the macrospin approx- 
imation adequately describes the system [34j. This is 
indeed the case: the free layer magnetization remains in 
a macrospin-like state for the studied range of I 28j . We 
show in Figs. [DandE^a) that the micromagnetic results 
produce power spectral densities (PSDs) and plane cross- 
ing angle distributions that are in good qualitative agree- 
ment with the macrospin simulations and experimental 
data. The computational demands of finite-temperature 
micromagnetic simulations preclude the use of a fitting 
procedure akin to that employed for macrospin simula- 
tions, thus we ran micromagnetic simulations with the 
parameters identified from the macrospin fits. The excep- 
tion is the magnitude of the ST polarization efficiency P, 
which was increased to 0.70 from 0.60 in the macrospin 
case. This discrepancy may stem from both a decrease of 
the magnetoresistive signal arising from non-uniformities 
in the free layer magnetization and additional dissipation 
of the energy supplied by ST drive via micromagnetic de- 
grees of freedom. Despite a small overall shift in the gen- 
eration frequency (« 0.3 GHz), we find that the PSDs 
given by micromagnetic simulations are in good agree- 
ment with the macrospin simulations and the experiment 
except at the lowest bias currents studied. Local mag- 
netization pinning at defects (discretization artifacts) or 
nonlinear damping [35J may be partially responsible for 
the discrepancy. 

While providing a direct means of extracting ST pa- 
rameters from the experiment, simulations of the stochas- 
tic LL equation give no insight into the mechanisms lead- 
ing to the observed asymmetries of the crossing angle 
distributions and spectral line sha pes . We turn to the 
magnetic Fokker-Planck equation [36|, which describes 
the deterministic drift and diffusion of the probability 
distribution of m on the unit sphere. Working in the en- 
ergy coordinate of the system instead of angles 9 and (p, 
one can derive 22|, [23j a Boltzmann-form probability per 
unit area of the sphere, p'(E), that the system possesses 
a particular energy, E: 



p'(E,I) = -exp(-P VE cS (E,I)). 



(5) 



Here V is the free layer volume, f3 = l/k^T, and Z is 
the partition function [11]. The ST-dependent effective 
energy E e g(E, I) is constructed by integrating the non- 
conservative torques acting on the magnetization along 



4 



— i ' 1 ' r 

Experiment 

Micromagnetic 

Macrospin 



1 ' 1 ' 1 ' — 

0.8 mA . ( a ) 




0.5n 
In-Plane Angle cp 



1.5n 



2.8 
2.6 
2.4 
2.2 
2 
1.8 



1(c) 


— ttt — i — 1 — ; 

o\mA \ - 








0.2mA\\ \ 

\\ 


.'j-''. 


i . i ,T- 



-0.15 -0.1 -0.05 
Energy E/(2nM;r) 




2.4 2.6 2.8 3 
Frequency (GHz) 



FIG. 3: (Color Online) (a) Measured PSDs delivered to a 
50Q load along with macrospin and micromagnitic simulation 
results (/ in steps of 0.1mA). (b) Effective energies E e g(<p) 
and corresponding crossing probability distributions p c (</>)■ 
(c) Frequencies (solid lines) and generated electrical powers 
(dotted lines) of the macrospin conservative orbits, (d) Com- 
parison of PSDs from macrospin simulations and the effective 
energy approach given by Eq.([7|). 



the conservative orbit of energy E. This approach re- 
lies on the assumption that the magnetization mainly 
evolves along conservative orbits, though it is induced by 
non-conservative and thermal torques to spread among 
the orbits on a time scale much longer than that of the 
oscillation period. 

For comparison with our measurements, we are inter- 
ested not in p'(E), but rather in the probability p c (f) 
that m crosses through the plane of the sample at angle 
tp. This quantity is given by 



Pc(<P) 



2Tr 1 M s pi(E( l p),I) 
Z u{E{fp)) 



dE(<p) 



dip 



(6) 



where uj(E) is the angular frequency and E(ip) is the 
in-plane cross-section of the conservative energy surface 
|28j . The distributions p c {tp) and energy surfaces E g(E) 
from which they are derived are shown in Fig. EJb). The 
Fokker-Planck approach necessarily predicts zero prob- 
ability at the minimum of E(ip) due to the vanishing 



density of states \dE((p)/dtp\ appearing in Eq. (J6j> . The 
crossing distributions are also plotted in Fig. [2] wherein 
we see excellent agreement with the macrospin results 
over the entire range of currents. We note that Eq. (J5j) 
and Eq. ([5]) can be easily inverted to reconstruct the ef- 
fective energy surface E C ^{E 1 I) from the measured plane 
crossing probability distributions p c (<p). 

We can now identify the cause of asymmetry and 
broadening of p c {f) observed at large values of /. While 
E e g((p) is approximately quadratic near the bottom of 
the well, it quickly crosses into a nearly linear regime 
with increasing ip. The number of available orbits be- 
tween E and E + dE becomes large in this region of 
E e s(ip) and m spends proportionally more time on these 
large amplitude trajectories. This causes the tail of p c (ip) 
to elongate. As the slope of this linear region decreases 
in response to increasing /, the distributions widen and 
distort commensurately. 

For currents near I c , flatness at the bottom of the ef- 
fective energy well E c g(ip) suggests a simple method for 
evaluating the STO's PSD. The relaxation of m towards 
its equilibrium orbit is driven by deterministic torques 
proportional to the slope of E e s(ip). Thus, near the bot- 
tom of E e g(ip), these restoring torques become small and 
the time evolution of the oscillation amplitude becomes 
dominated by thermal diffusion. If the time scale for 
thermal diffusion of the amplitude is long compared to 
the period of oscillations, the PSD can be approximated 
by a superposition of the power generated by each orbit 
weighted by the probability of finding m on this orbit 
given by the simple Boltzmann-like expression of Eq. ([5]) . 
We proceed with this simple method for determining the 
PSD, and thereby verify that this diffusion-dominated 
limit is applicable for our system. We first compute the 
frequencies w{E) and average electrical powers P(E) (de- 
livered to a 5051 load) corresponding to each conservative 
orbit, both of which are plotted in Fig. |3jc). For a single- 
valued function E(ui), the PSD is given by the electrical 
power of individual orbit P{E) weighted by the proba- 
bility of finding the system at this orbit: 



S(cu) 



2irjM s p'(E(uj),I)P(E(u])) 



(J) 



For a multi-valued function E(uS), such as that shown in 
Fig. G3c), a sum over all branches of E(uj) should be 
added to the right hand side of Eq. ([7]). As shown in 
Fig. Eld), the PSDs calculated with Eq.© possess pow- 
ers, line widths, and line shape asymmetries that are in 
excellent agreement with the macrospin simulations and 
experimental PSDs for / close to I c . Far from I c , the 
agreement between PSD given by the macrospin simu- 
lations and by Eq.([7]) decreases because rapid amplitude 
relaxation and an associated increase in phase noise (aris- 
ing from nonlinear coupling between oscillation ampli- 
tude and phase 15j ) belie the simplicity of our ensemble 



average. 
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With Eq. (JT]) we can now see that the PSD broad- 
ening and asymmetry near I c has the same origin as the 
asymmetry of the crossing angle distributions p c {ip). The 
occupation of larger amplitude orbits increases substan- 
tially as Eeg(tp) develops a large linear region. As seen 
in Fig. [3J these orbits posses lower frequencies due to the 
frequency redshift with increasing amplitude. Thus the 
PSDs are seen to spread asymmetrically towards lower 
frequencies. In principle, non-monotonicities in lo(E) 
such as those seen in Fig. [3jb) can create, by means of 
an increased density of states, an accumulation of power 
near local extrema of ou(E). Closely spaced features of 
this sort could even give the illusion of multiple modes. 
Since our samples are prone to dielectric breakdown at 
higher current densities, we are unable to experimentally 
access this regime. 

In conclusion, time domain measurements of STO volt- 
age allow us to rapidly map statistical ensembles of STO 
magnetization trajectories and thereby determine ST- 
dcpcndcnt Fokker-Planck effective energy of the STO. 
We use the ST-dependent effective energy approach to 
derive a simple expression for the STO power spectral 
density valid at / w I c , which is in excellent agreement 
with our experimental data. Our technique also provides 
a simple and accurate measurement of the field-like com- 
ponent of ST across a wide range of voltage bias values. 
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No. HR0011-10-C-0153 and by NSF Grants No. DMR- 
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RECONSTRUCTING THE REAL-SPACE 
DYNAMICS 

Accurate determination of the spatial orientation of 
the free layer magnetization m is predicated upon care- 
ful characterization of a sample's voltage-dependent tun- 
neling magnetoresistance (TMR). To this end, the max- 
imum and minimum observed resistances Rap(V) and 
Rp (V) are extracted as a function of bias from hysteresis 
loops taken in the vicinity of the sample's easy axis. By 
measuring within a ±5° range of angles we are able to 
account for potential field misalignments. Representative 
hysteresis loops are shown for both voltage polarities in 
Figs. (Ha) andQJb), wherein the influence of magnetiza- 
tion dynamics on the average resistance of the sample is 
readily visible (positive current pulls the magnetization 
towards the parallel state in this geometry). We have 
confirmed that the powers of these magnetization oscil- 
lations are negligible at high fields (approaching 1.5 kOe), 
and therefore conclude that the opening angles of the tra- 
jectories are too small to interfere with measurements of 
Rp(V) and R A p{V). 

The bounding resistance values, as well as the inverse 
of the average conductance Gq 1 (corresponding to an an- 
gle ip — tt/2 between m and p) are plotted in Fig. [TTc) 
as a function of bias. These traces provide all of the nec- 
essary information for mapping the sample's resistance 
state back to its orientation by means of Eq. 3 of the 
main text. 



Obtaining the time-dependent sample resistance still 
requires some care. We must determine the voltage os- 
cillations across the sample AV(t) corresponding to the 
signal V(t) measured at the oscilloscope. This requires 
characterization of our microwave circuit as detailed in 
Fig. [2] After accounting for active amplification and 
passive circuit attenuation, the desired signal is given by 



AV(t) = V(t) 



5on + R e 



(R) 



5on 



(i) 



where i? ox = 210 is the resistance contribution from 
sources extrinsic to the sample itself (from contacts and 
probes), and (R) is the time averaged MTJ resistance 
(less the extrinsic contribution) that we measure using 
the sourcemeter that supplies direct current to the de- 
vice. The voltage oscillations measured by the oscillo- 
scope may then be mapped back to the angular separa- 
tion of the free layer and polarizer. As described in the 
main text, we can only pinpoint the location of the mag- 
netization as it crosses the plane of the sample, which 
corresponds to extrema in R(t). It remains, then, to 
develop a robust peak-finding procedure in order to cor- 
rectly identify the extrema from R(t). 

The primary culprits of reconstruction errors are er- 
roneous extrema originating from any source of high- 
frequency low-amplitude noise, such as thermal fluctu- 
ations of the magnetization orientation or an extrinsic 
source such as electronics noise. Such anomalies are eas- 
ily eliminated by the use of a Fourier smoothing proce- 
dure that eliminates high frequency features in the volt- 
age signal. We low-pass filter all frequency components 
greater than 1.5 times the fundamental frequency of the 
oscillator, beyond which the experimentally measured 
power spectral density is nearly zero. We apply a mag- 
netic field of sufficient strength and appropriate direction 
to intentionally constrain the ensemble of orbits away 
from the easy axis. Otherwise a second harmonic would 
be present in the voltage signal, leading to a more compli- 
cated crossing angle measurement procedure. Also, since 
low-frequency noise can smear out the reconstructed an- 
gular distributions, we only admit frequencies above 300 
MHz, below which point the noise figure of our ampli- 
fier rapidly increases from its base value of 1.3dB. We do 
not observe any large low-frequency features in the ex- 
perimentally obtained spectra of the STO, and therefore 
do not expect that this high-pass filtering will impact 
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FIG. 1: (Color Online) Resistance as a function of field at various values of (a) positive and (b) negative bias voltage. The 
increment between adjacent traces is 50mV. (c) The corresponding Rap(V), Rp(V), and l/Go(V) as extracted from (a) and 
(b). The extrinsic resistance contribution i? cxt is included in (a-c), though it is subtracted for the actual analysis. 




FIG. 2: (Color Online) The microwave circuit used for time 
domain analysis. The same instrument is used to supply cur- 
rent I and measure the time averaged resistance across the 
DC branch of the circuit. The amplified microwave voltage 
generated by the sample, AV(t), results in the voltage V(t) 
across the characteristic 50O real time oscilloscope. 

the accuracy of the method. Examples of both raw and 
band-pass filtered time-domain data are shown in Fig. [3] 
(a, b) along with an example of how the low-pass filtering 
affects the crossing angle distributions (Fig. [3jc)). 

A severe micromagnetic curvature of the free layer 
would preclude use of our method, since the projection 
m • p would no longer provide a complete representa- 
tion of the free layer's state. Later on in this supplement 
we present the results of micromagnetic simulations that 
suggest there is minimal distortion of the magnetization 
profile from the macrospin state. 

A broadening of the reconstructed peaks will inevitably 
result both from small micromagnetic distortions and 
from the noise at the oscilloscope. We therefore expect 



that the angular distributions that we obtain are actually 
convolutions of the intrinsic distributions with an approx- 
imately Gaussian kernel. Deconvolution is unavailing in 
this particular context, so we stress that the distribu- 
tions we plot are slightly smeared versions of the true 
ones. This effect is most pronounced for small applied 
currents since the relative contribution of the electronics 
noise (« 5 mV rms ) is largest in that case. At most we 
expect a 3-4° blurring of the distributions, which does 
not alter any of their qualitative features. 



MACROSPIN SIMULATIONS 

Macrospin simulations allow us to make quick com- 
parisons of our experimental results to theoretical expec- 
tations. We adopt the stochastic Landau-Lifshitz (LL) 
form of the magnetization dynamics 

dm = [v(m, t) — v 2 m\ dt — m x (z^dW + avm x dW) , 

(2) 

interpreted in the Ito sense [2j . The deterministic part of 
the dynamics is described by 

v(m, t) = m x (h off + am x (h cff + /3 st p/a)) , (3) 

where a is the Gilbert damping parameter, v is the mag- 
nitude of thermal fluctuations as given by the fluctuation- 
dissipation theorem, and dW is the isotropic Wiener pro- 
cess that is the generator of Brownian motion. All fields 
and magnetic moments are normalized by the saturation 
magnetization M s . In Eqs. @ and ([3]) only, the units 
of time are are multiplied by ^M s (where 7 is the gyro- 
magnetic ratio) to yield a dimensionless quantity. The 
spin-torque coefficient /3 s t depends on the mutual angle 
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FIG. 3: (Color Online) Two representative segments (a,b) of 
raw and band-pass filtered ( "Fourier smoothed" ) time-domain 
data. The horizontal dotted line represents the time-averaged 
resistance (R). Local extrema are seen to be effectively re- 
moved by low-pass filtering. The effects of high-pass filtering 
are shown in (c) where the reconstructed crossing distribu- 
tions are shown with and without such a filter. All data are 
shown for I = 0.7mA, H = 500Oe, and ip = 120°. 



of m and p and is given by [3] 



P 



2eM s 2 V Vl + ^ 2 m-p 



(4) 



hat = 



2eM s 2 V Vl + ^ 2 m-p 



(5) 



is taken to be a constituent of the effective field h c ff. 
Here, a^t, represents the dimensionless strength of the 
field-like torque. We integrate Eq. ([2]) over many real- 
izations of the Wiener process simultaneously on a GPU 
using an Euler algorithm with appropriately small time 
increment At = 30. x 1CP 15 s. The probability distribu- 
tions calculated from these newly alacritous simulation 
results can be fed directly into a least-squares optimiza- 
tion algorithm. The simulation results are interpolated 
onto the same set of cp points as the experimental distri- 
butions, and the residuals between the two are minimized 
with respect to a s tt, Oflt, and the first-order perpendicu- 
lar magnetic anisotropy field H p q. The plots of the FLT 
strength in the main text is given, instead, in units of 
field by utilizing the whole prefactor of Eq. ([5]). 

The spectral properties of these dynamics are quanti- 
ties of great interest. In order to compare to experimen- 
tally observed power spectral densities, we must calculate 
the actual voltage signal that the simulated device would 
deliver, ceteris parabis, to a 50il load in the experimental 
circuit. This procedure is somewhat complicated by the 
nature of the measurement: in the presence of a constant 
current / across the device, changes in the orientation 
of the free layer alter the voltage across the junction, 
which in turn changes the observed resistance through 
the bias-dependence of TMR. We must therefore derive 
an expression for the voltage as a function of both / and 
m. Starting from the expression for the conductance as 
a function of angle and voltage, 



G(V,9) = y=G (V) 



1 



AG(U)m-p, 



(G) 



and approximating the functions Go(V) and G{V) as 

Gq(V) = 9 o+9iV + g 2 V 2 (7a) 
AG(V) = h + h 1 V + h 2 V 2 , (7b) 

we solve this system of equations for V and expand in 
powers of cos(#) and I to whichever order is necessary 
for a reasonable facsimile to the full solution. Using the 
final approximation of V(I, m) we calculate the voltage 
signal at N individual points Vi separated by time At. 
These points are then multiplied by a Hann window in 
order to mitigate spectral artifacts caused by the finite 
signal lengths. The RMS power delivered to a 50J7 load 
at spectral component k is then 



where P is the polarization efficiency, / is the current, 
and V is the free layer's volume. The dimensionless factor 
dst gives the strength of the in-plane spin torque, and is 
one of the fitting parameters used when comparing the 
simulation results to experimental data. The field-like 
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N 



50fl + R c 



(R) 
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where T[xi] represents the (un-normalized, single-sided) 
discrete Fourier transform, Wi are the sampled points of 
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the Hann windowing function, and the factor of 2 in 
the numerator accounts for the power lost in the first 
harmonic from the use of the windowing function. The 
power spectra are calculated separately for each realiza- 
tion of the thermal field, and are then averaged to obtain 
the final spectrum. 



EXTRACTING THE FIELD-LIKE TORQUE 

Extraction the field-like torque from crossing distri- 
butions involves many fewer complications than deriv- 
ing the same results from ST ferromagnetic resonance 
(STFMR) methods. We assume a Stoner-Wohlfarth 
(SW) energy density of the form 

E{tp) = 2vr Ml (N xx - N yy ) cos 2 <p 

- (H x + H & t)M s cos p - H y M s sin ip, (9) 

where H x and H y are the in-plane components of the 
external field (including the dipolar contribution H^ ip ), 
N xx and Nyy are the in-plane components of the demag- 
netizing tensor, and Hfn is the FLT contribution we seek 
to identify. We find the equilibrium azimuthal angle ipo 
that minimizes Eq. ^ with iJfl t = 0, then we express 
the magnitude of the FLT 

H&t = -H x + H y cot(<y9 + Sip) 

+4ttM s (N xx - iVyy) cos^o + S<p) (10) 

in terms of the FLT-induced deflection from equilibrium 
Sip. By comparing the displacements, rather than abso- 
lute angles, we can partially account for errors (discussed 
below) induced a slight micromagnetic curvature of the 
system. It would be natural to compare the results fur- 
nished by this method to those from STFMR, though we 
were unable to do so for our sample given the presence 
of large backgrounds in our finite-bias resonance spectra 
(a typical complication of such measurements). 



MICROMAGNETIC SIMULATIONS 

A more realistic representation of the experimental sys- 
tem is obtained in the micromagnetic approximation. To 
this end we modify the OOMMF finite-difference micro- 
magnetic framework to perform finite temperature simu- 
lations of the LL equation including spin-torque in a form 
appropriate for MTJs The computational burden 
of these simulations is great: Gaussian random deviates 
must now be generated for each site in the discretized 
system, and the time step of the simulations must be 
drastically reduced (from the zero temperature case) in 
order that the timescale of thermal fluctuation remains 
much smaller than the timescale of the magnetization 
dynamics. 
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FIG. 4: (Color Online) Snapshot of the micromagnetic config- 
uration of the sample at I — 0.8mA (the largest current value 
used in this study), wherein color encodes the out-of- plane 
component of magnetization. 



For these simulations we set the exchange stiffness A 
to 2.0 x 10~ 6 emu/cm, while the remaining simulation 
parameters are identical to those used in the macrospin 
approximation. We deliberately omit a realistic magneto- 
static field from the polarizer in favor of a simple uniform 
field for the sake of a more straightforward comparison 
to the macrospin results. The aim of this work is not to 
reproduce the exact local environment of the experimen- 
tal system, but rather to gauge the agreement of these 
various numerical and analytical methods. 

During the course of these simulations the average 
magnitude of the magnetization unit vector (as defined 
by the norm of the vector average over cells) spends 
approximately 50% of its time above 0.95 and approx- 
imately 90% of its time above 0.90. This is a good in- 
dication that the micromagnetic curvature of the system 
has a minimal impact on our reconstruction process. We 
show in Fig. [4j for reference, a high-curvature snapshot 
of the magnetization along a large amplitude trajectory. 

The micromagnetic simulations resolve an apparent 
inconsistency between the experiment results and the 
macrospin simulations: the equilibrium angles predicted 
within the macrospin approximation (i.e. by the Stoner- 
Wolfarth model) do not match those observed in experi- 
ment; they differ by a factor of 0.037T. This discrepancy 
is due to a slight micromagnetic curvature of the sys- 
tem, though we emphasize that this curvature found in 
static equilibrium is even smaller than dynamic curvature 
observed during large-amplitude motion of the magneti- 
zation pictured in Fig. |4l 



THE EFFECTIVE ENERGY APPROACH 

The aforementioned simulation techniques do little to 
further our understanding of the system beyond what is 
already apparent from measurements. We seek to under- 
stand, primarily, the origin of asymmetry and broaden- 
ing in both the ip distributions and spectral lines. An 
analytical approach is therefore sought. Since we are di- 
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FIG. 5: (Color Online) (a) Logarithmic and (b) linear scale plots of the energy distributions according to Eq. (|16[) . The 
distributions are seen to diverge slowly from the Boltzmann result at / = 0. The minimum available energies at successive 
currents are observed to shift. This is due to modifications of the energy surface by FLT and voltage induced anisotropy, 
which together decrease the lowest attainable energy value with increasing I. One can also observe small peaks in the high-bias 
distributions just above the energy minima. These are indications that the system has reached the critical current and favors 
an orbit of non-zero radius corresponding to an energy above that at the bottom of the well. 



rectly measuring probability distributions, it behooves us 
to leave behind the stochastic differential equation Eq. 
@ in favor of the equivalent Fokker-Planck equation. 
Following the original derivation of the magnetic Fokker- 
Planck equation [5J , it has since been recast in the effec- 
tive energy coordinate of the system [l|, Q . The central 
idea of this approach is to assume that non-equilibrium 
and thermal torques cause the system to drift slowly be- 
tween conservative orbits of the system (which are by 
definition calculated in the absence of such drives). For 
a magnetic system with a single potential well, the effec- 
tive energy associated with a conservative orbit of energy 
density E is 

E cS (E) - £ cff (£o) = {E-E )-I r){E')dE', (11) 

J Eq 

where Eq is an arbitrary reference energy which we 
choose to be the maximum of the in-plane energy E(tp). 
The factor Ii](E) is the ratio of the works done by ST 
and damping over the conservative orbit of energy E. 
The current / has been factored out in order to show 
the explicit current dependence. The functional form of 
r](E) can be found from the constituent torques of Eq. 
([2]). Aside from energy contributions of the demagne- 
tizing and external fields, we include those due to per- 
pendicular anisotropy and field-like torque (FLT) , which 
we treat as an effective field contribution. The energy 



density from perpendicular anisotropy is given by 

E an = -^sm 2 9 ~ K 2 siu i 9 1 (12) 

where K\ and K 2 are the first and second order 
anisotropy energies: 

Ki = ^(H p0 + AH p0 )M s (13) 

K 2 = - A H pl M s . (14) 

In addition to the first and second order anisotropy fields, 
H p Q and H p i , we include the voltage-controlled magnetic 
anisotropy (VCMA) contribution AH p o which we take to 
be 600 Oe/V based on experiments conducted on similar 
samples [7[. The field-like torque's energy contribution 
is, meanwhile, given by [l[ 

^ = -|^ln(l + P 2 m .p). (15) 

When including either of these voltage and current- 
dependent terms we must recalculate the set of conserva- 
tive orbits at each current step since they modify the ef- 
fective field. We employ our aforementioned GPU-based 
macrospin simulation code to simultaneously calculate 
the entire set of conservative orbits. The probability per 
unit area of the unit sphere of finding the magnetization 
at energy E is given in the convenient Boltzmann form: 

p'(E, I) = 1 exp(-/3 VE cS (E, I)). (16) 
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where j3 is the inverse of the thermal energy (fe^T) -1 , 
and the current has been added as an argument of the 
effective energy surface since it must be recalculated for 
every choice of I. The factor 



iM, 



f B l t{E') P '{E',I)&E' 



(17) 



is a "non-equilibrium partition function" resulting from 
an integral over all energy states (see Eq. 10.203 and 
related discussion of Ref. [l|). Here, t(E) are the pe- 
riods corresponding to conservative orbits of energy E. 
The p'(E) distributions are shown for several values of 
the current in Fig. [5j where one can readily observe the 
transition from Boltzmann to non-Boltzmann behavior 
with increasing I. We are interested in the probability 
that the magnetization crosses through the plane of the 
sample at angle ip. This quantity is simply related to 
p'{E), and is found through the "conservation of proba- 
bility" to be 



p c (ip) = ±p'(E(ip),I)A(E(ip)) 



dE{ip) 



dip 



(18) 



In this expression, the density of states accounts for the 
number of orbits between ip and dip, while E(ip) de- 
notes the in-plane cross-section of the energy surface 
with 9 = 7r/2. The crossing probability p c {ip) always 
goes to zero at the center of the well since the density 
of states \dE(ip)/dip\ vanishes at that point. The fac- 
tor A(E)dE = ^M s r(E)dE gives the surface area of the 
unit sphere enclosed by the band of orbits between E 
and E + dE, and is proportional to the period of the or- 
bits t{E) in that vicinity (as shown in the supplemental 
material of Ref. @) . The area factor, a minor correction 
for our system, is necessary since p'{E) is specified per 
unit area of the unit sphere. We may rewrite the crossing 
probability in a more intuitive form as 



Pc{ip) 



2tt 1 M sP '{E( V ),I) 
Z u{E(ipj) 



dE(ip) 



dip 



(19) 



where uj(E) is the angular frequency. Thus, we need only 
record uj(E) and P(E) (the power dependence) in order 
to furnish predictions for the spectral line shapes of the 
STO. This is an important result, as such quantities are 
readily obtained from simulations. 

In the approximation of this model (as described in 
the main text), the line shape is given by the sum of 
the powers of individual orbits weighted by their relative 
occupation probabilities: 



S(u) = p'(E(u),I)A(E( U ))P(E(u)). 



(20) 



We can make a transformation from A(E) to 
2ir"/M s /ijj(E) to yield the expression in the main text: 



ui 



(21) 



The application of this method to the experimental sys- 
tem described in the main text has proved very successful 
near threshold. It remains that we should compare this 
method to a well-established formalism that is widely 
used to describe STO dynamics. 



COMPARISON TO NONLINEAR 
AUTO-OSCILLATOR THEORY 




FIG. 6: (Color Online) Comparison of the linewidths pre- 
dicted by the full KTS auto-oscillator formalism (solid blue 
lines) along with their asymptotes (dotted blue lines), and 
the results from applying Eq. (|21[1 (solid red lines). N is a 
measure of the system's nonlinearity, negative values of which 
indicate a redshift with increasing current /. The abscissa is 
scaled by the critical current 7 C , while the vertical axis is given 
in terms of the half-width at half-maximum Aoj and the linear 
relaxation rate Tq- Data from the full KTS formalism was 
taken from Fig. 3 of Ref. @- 

The nonlinear auto-oscillator theory developed by 
Kim, Tiberkevich, and Slavin (KTS) has been quite suc- 
cessful at predicting spectral properties of STOs §£1. 
The caveat of the KTS approach is that while it predicts 
simple asymptotic forms for the sub- and super-threshold 
linewidths, results near threshold can only be calculated 
by the involved process of finding the eigenfunctions and 
eigenvalues of the Fokker-Planck equation for the oscil- 
lator's power. Our approach, meanwhile, performs best 
precisely in the realm where the KTS theory becomes 
computationally burdensome. 

We plot in Fig. |6]the KTS results for the generation 
linewidth of the magnetic system described in Ref. H|. 
We then use the stationary probability distributions (as 



per Eq. (83b) of Ref. [bj), dimensionless powers, and 



frequencies calculated within the KTS framework (ignor- 
ing the explicit character of the underlying conservative 
orbits) as inputs to Eq. (|2"Tj) in order to judge the agree- 
ment with our ensemble-averaged method. There is a 
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slight difference in the plotted quantities: while the KTS 
linewidths are taken from a single Lorentzian fit to an 
asymmetric power spectral density, our own results are 
simply the half-width at half-maximum taken from the 
distributions calculated using Eq. (P?Tj) . Our method fa- 
vorably reproduces the spectral line widths near thresh- 
old, confirming its generality. 

Combining our expression with the asymptotes of the 
KTS formalism, one can now easily calculate the ex- 
pected linewidths of a STO across the entire domain of 
currents from sub- to super-threshold. 
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